{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using VMLS\n",
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 8\n",
    "# Linear equations\n",
    "### 8.1 Linear and affine functions\n",
    "**Matrix-vector product function.** Let’s define an instance of the matrix-vector\n",
    "product function, and then numerically check that superposition holds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3 Array{Float64,2}:\n",
       " -0.1   2.8  -1.6\n",
       "  2.3  -0.6  -3.6"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [-0.1 2.8 -1.6; 2.3 -0.6 -3.6] # Define 2x3 matrix A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "f (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(x) = A*x # Define matrix-vector product function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([9.47, 16.75], [9.47, 16.75])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let’s check superposition\n",
    "x = [1, 2, 3]; y= [-3, -1, 2];\n",
    "alpha = 0.5; beta = -1.6;\n",
    "lhs = f(alpha*x+beta*y)\n",
    "rhs = alpha*f(x)+beta*f(y)\n",
    "lhs,rhs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.7763568394002505e-15, [2.8, -0.6], [2.8, -0.6])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(lhs-rhs), f([0,1,0]), A[:,2] # Should be second column of A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**De-meaning matrix.** Let’s create a de-meaning matrix, and check that it works\n",
    "on a vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([-0.966667, 1.13333, -0.166667], [-0.966667, 1.13333, -0.166667])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "de_mean(n) = eye(n) .- 1/n; # De-meaning matrix\n",
    " x = [0.2, 2.3, 1.0];\n",
    " ex1 = de_mean(length(x))*x # De-mean using matrix multiplication\n",
    " ex2 = x .- avg(x) # De-mean by subtracting mean\n",
    "ex1, ex2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Examples of functions that are not linear.** The componentwise absolute value\n",
    "and the sort function are examples of nonlinear functions. These functions are\n",
    "easily computed by `abs` and `sort`. By default, the `sort` function sorts in increasing\n",
    "order, but this can be changed by adding an optional keyword argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "f (generic function with 1 method)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(x) = abs.(x) # componentwise absolute value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Int64,1}:\n",
       " 1\n",
       " 2"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = [1, 0]; y = [0, 1]; alpha = -1; beta = 2;\n",
    "f(alpha*x + beta*y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Int64,1}:\n",
       " -1\n",
       "  2"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha*f(x) + beta*f(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Int64,1}:\n",
       "  2\n",
       " -1"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(x) = sort(x, rev = true) # sort in decreasing order\n",
    "f(alpha*x + beta*y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Int64,1}:\n",
       " 1\n",
       " 0"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha*f(x) + beta*f(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.2 Linear function models\n",
    "**Price elasticity of demand.** Let’s use a price elasticity of demand matrix to pre-\n",
    "dict the demand for three products when the prices are changed a bit. Using this\n",
    "we can predict the change in total profit, given the manufacturing costs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "77.51999999999998"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = [10, 20, 15]; # Current prices\n",
    "d = [5.6, 1.5, 8.6]; # Current demand (say in thousands)\n",
    "c = [6.5, 11.2, 9.8]; # Cost to manufacture\n",
    "profit = (p-c)'*d # Current total profit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " -0.3   0.1   -0.1 \n",
       "  0.1  -0.5    0.05\n",
       " -0.1   0.05  -0.4 "
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Demand elasticity matrix\n",
    "E = [-0.3 0.1 -0.1; 0.1 -0.5 0.05 ; -0.1 0.05 -0.4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " -0.1                \n",
       "  0.05               \n",
       " -0.06666666666666667"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_new = [9, 21, 14]; # Proposed new prices\n",
    "delta_p = (p_new-p)./p # Fractional change in prices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "  0.04166666666666667\n",
       " -0.03833333333333334\n",
       "  0.03916666666666667"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta_d = E*delta_p \n",
    "# Predicted fractional change in demand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 5.833333333333333\n",
       " 1.4425           \n",
       " 8.936833333333333"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_new = d .* (1 .+ delta_d) \n",
    "# Predicted new demand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(66.25453333333333, 77.51999999999998)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "profit_new = (p_new-c)'*d_new \n",
    "# Predicted new profit\n",
    "\n",
    "profit_new, profit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we trust the linear demand elasticity model, we should not make these price\n",
    "changes.\n",
    "\n",
    "**Taylor approximation.** Consider the nonlinear function $f : R^2 → R^2$ given by\n",
    "$$\n",
    "f(x) =\n",
    "\\begin{bmatrix}\n",
    "‖x − a‖\\\\\n",
    "‖x − b‖\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix} \n",
    "\\sqrt{(x_1 − a_1)^2 + (x_2 − a_2)^2}\\\\\n",
    "\\sqrt{(x_1 − b_1)^2 + (x_2 − b_2)^2}\n",
    "\\end{bmatrix}\n",
    ".\n",
    "$$\n",
    "\n",
    "The two components of $f$ give the distance of $x$ to the points $a$ and $b$. The function is differentiable, except when $x = a$ or $x = b$. Its derivative or $Jacobian$ matrix is given by\n",
    "\n",
    "$$\n",
    "Df(z) =\n",
    "\\begin{bmatrix}\n",
    "\\frac{∂f_1}{∂x_1}(z) & \\frac{∂f_1}{∂x_2}(z) \\\\\n",
    "\\frac{∂f_2}{∂x_1}(z) & \\frac{∂f_2}{∂x_2}(z)\n",
    "\\end{bmatrix}\n",
    " = \n",
    "\\begin{bmatrix}\n",
    "\\frac{‖z_1 − a_1‖}{‖z − a‖} & \\frac{‖z_2 − a_2‖}{‖z − a‖} \\\\\n",
    "\\frac{‖z_1 − b_1‖}{‖z − b‖} & \\frac{‖z_2 − b_2‖}{‖z − b‖} \n",
    "\\end{bmatrix}\n",
    ".\n",
    "$$\n",
    "Let’s form the $Taylor$ approximation of $f$ for some specific values of $a$, $b$, and $z$,\n",
    "and then check it against the true value of $f$ at a few points near $z$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 0.9055385138137417\n",
       " 1.2727922061357855"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(x) = [ norm(x-a), norm(x-b) ];\n",
    "Df(z) = [ (z-a)' / norm(z-a) ; (z-b)' / norm(z-b) ];\n",
    "f_hat(x) = f(z) + Df(z)*(x-z);\n",
    "a = [1, 0]; b = [1, 1]; z = [0, 0];\n",
    "f([0.1, 0.1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 0.9               \n",
       " 1.2727922061357857"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_hat([0.1, 0.1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 0.7071067811865476\n",
       " 0.7071067811865476"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f([0.5, 0.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 0.5               \n",
       " 0.7071067811865477"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_hat([0.5, 0.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regression model. We revisit the regression model for the house sales data in\n",
    "Section [2.3](https://web.stanford.edu/~boyd/vmls/vmls.pdf#section*.192). The model is\n",
    "\n",
    "$$\n",
    "ŷ = x^Tβ + v = β_1x_1 + β_2x_2 + v,\n",
    "$$\n",
    "\n",
    "where $ŷ$ is the predicted house sale price, $x_1$ is the house area in $1000$ square feet, and $x_2$ is the number of bedrooms. \n",
    "\n",
    "In the following code we construct the $2 × 774$ data matrix $X$ and vector of outcomes $y^d$, for the $N = 774$ examples in the data set. We then calculate the regression model predictions $ŷ^d$, the prediction errors $r^d$, and the $RMS$ prediction error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(774, (2, 774))"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# parameters in regression model\n",
    "beta = [148.73, -18.85]; v = 54.40;\n",
    "D = house_sales_data();\n",
    "yd = D[\"price\"]; # vector of outcomes\n",
    "N = length(yd)\n",
    "X = [ D[\"area\"] D[\"beds\"] ]';\n",
    "N, size(X) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "74.84571862623022"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ydhat = X'*beta .+ v; # vector of predicted outcomes\n",
    "rd = yd - ydhat; # vector of predicted errors\n",
    "rms(rd) # RMS prediction error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "112.78216159756509"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compare with standard deviation of prices\n",
    "stdev(yd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3 Systems of linear equations\n",
    "**Balancing chemical reactions.** We verify the linear balancing equations on page [155](https://web.stanford.edu/~boyd/vmls/vmls.pdf#section*.192)\n",
    "of VMLS, for the simple example of electrolysis of water."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Int64,1}:\n",
       " 0\n",
       " 0"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = [2 ; 1]\n",
    "P = [2 0 ; 0 2]\n",
    "# Check balancing coefficients [2,2,1]\n",
    "coeff = [2,2,1];\n",
    "[R -P]*coeff"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.1.1",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
